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Abstract 



We study binary state contagion dynamics on a social network where nodes act in response to 
the average state of their neighborhood. We model the competing tendencies of imitation and non- 
conformity by incorporating an off-threshold into standard threshold models of behavior. In this way, 
we attempt to capture important aspects of fashions and general societal trends. Allowing varying 
amounts of stochasticity in both the network and node responses, we find different outcomes in the 
random and deterministic versions of the model. In the limit of a large, dense network, however, we 
show that these dynamics coincide. The dynamical behavior of the system ranges from steady state 
to chaotic depending on network connectivity and update synchronicity. We construct a mean field 
theory for general random networks. In the undirected case, the mean field theory predicts that 
the dynamics on the network are a smoothed version of the average node response dynamics. Wc 
compare our theory to extensive simulations on Poisson random graphs with node responses that 
average to the chaotic tent map. 
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1 Introduction 



Almost universally, people enjoy observing, speculating, and arguing about their fellow humans' 
behavior, including what they are wearing, the music they listen to, and however else they express 
themselves. Fashions, trends, and fads intrigue us, whether past, present, or future. Contemporary 
examples of prevalent trends include skinny jeans, fixed- wheel bicycles, and products made by Apple 
Inc. In the 1960's, shag carpets and floral wallpaper would take their place. Trends are also 
prevalent in language: "hipster" was popular in the 1940's and evolved into "hippie" , but it has 
been reappropriated today with different connotations; "groovy" was once used seriously, and now 
comes out tongue-in-cheek. 

The dynamics of these cultural phenomena are fascinating and complex. They reflect numerous 
factors such as the political climate, social norms, technology, marketing, and history. These phe- 
nomena are also influenced by essentially random events. Usually, we can explain the emergence of a 
fad after the fact, but it is extremely difficult to predict a priori whether some behavior will become 
popular. Furthermore, those behaviors which are adopted by a large fraction of the population can 
lose their excitement and die out forever or later recur unpredictably. 

In this work we describe a mathematical model of the rise and fall of trends. In particular, we 
model a social contagion process where people are influenced by the behavior of their friends. The 
agents in the model act according to simple competing tendencies of imitation and non-conformity. 
One can argue that t hese two ingredients are essential to all trends; indeed, Simmel, in his classic es- 



say "Fashion" (| 1957T ) , believed that these are the essential forces behind the creation and destruction 
of fashions. 

Our model is not meant to be quantitative, except perhaps in carefully designed experiments, 
but it captures the features with which we are familiar: some trends take off and some do not, 
and some trends are stable while oth ers v ary wildly th r ough time. Our model is closely related to 



the seminal work of 



Schellind dl97lh and 



also c onnected to theor ies of percolatio n ( Stauffer and 



20031 ). and magnetism (jNewman . 



2003 



Granovetter ( 1978 ). Being a mathematical m odel, it is 



Aldana et al 



Ah aronv 



20031 ). 



19941 ) , disease spreading (jNewman 



We focus on the derivation and analysis of dynamical master equations that describe the expected 
evolution of the system state. Through these equations we are able to predict the behavior of the 
social contagion process. In some cases, this can be done by hand, but most of the time we resort 
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to numerical methods for their iteration or solution. 

This thesis is structured as follows. In Section [5J we introduce the reader to important back- 
ground material relevant to the on-off threshold model. In Section [31 we define the model and its 
deterministic and stochastic variants. In Section 01 we provide an analysis of the model when the 
underlying network is fixed. Section [5] develops a mean field theory of the model in the most general 
kind of random graphs. In Section [BJ we consider the model on Poisson random graphs with a spe- 
cific kind of response function. Our analysis is then applied to this specific case, and we compare the 
results of simulations and theory. Finally, Section [7] presents conclusions and directions for further 
research. 



2 Background 

2.1 Basics of graphs and networks 

When modeling any dynamical system of many interacting particles or agents, we are often forced 
to start with a simplified description of their interactions. In a solid, for example, atoms are situated 
on some sort of lattice and assumed to interact only with their nearest neighbors. However, in many 
cases the interactions we aim to model do not have the periodic structure of a lattice. Graphs, 
which are just a set of points connected by lines, are a more general mathematical structure. We 
denote a graph <S as an ordered pair ^ = (V, E) where V is the set of vertices (also called nodes) 
and E C V x V the set of edges (also called links). Here, E is a set of ordered pairs, and we denote 
an edge from vertices i to j as ij 6 E. The te rms graph and network will be used interchangeably, 



but they have acquired different connotations. 



Hacked ( 2011 



jl.2.1) puts this rather nicely: 



To model a complex system as a graph is to filter out the functional details of each 
of its components, and the idiosyncrasies of their interactions with each other, and to 
focus instead on the underlying structure (topology) as an inert mathematical construct. 
Although this technique is central also to network theory, the word network, in contrast, 
usually carries with it connotations of the context in which the overarching system exists, 
particularly when that system displays any sort of nonlinear dynamics. For example, 
when investigating the spread of infectious disease on a human sexual contact network 
it makes sense to consider the relevant sociological parameters as well as the abstract 
topology, and it is in such settings that the interdisciplinary aspect that distinguishes 
network theory comes to the fore. 

In models of human behavior, interactions can be considered to occur on a social network. Each 
person is connected to those people they interact with. Interactions that constitute a connection 
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Figure 1: A graph consisting of a connected component of 5 nodes and one isolated node. The 
majority of the edges are directed. There is one undirected edge, shown without arrows. There is 
one self-loop which should be interpreted as directed. 



between two people, A and B, may be defined in many ways. For example: 
Case 1. A contacts B. This could mean: 

(a) A sends B an email. 

(b) A attends the same concert as B. 

Case 2. A is B's superior in a hierarchy. 

Case 3. A and B both belong to the same group. 

These examples illustrate different types of edges that can be included in a graph. In Case Hal 
sending an email message is a one-way communication, so it forms a directed edge. We represent 
directed edges with arrows in a drawing of the corresponding graph. In Case llbl however, attending 
a concert is a symmetric relation and there is no reason to give the edge direction. Such edges 
are called undirected or bidirectional. Sec Figure [T] for an example of a small graph containing 
multiple types of edges. A simple graph, by definition, contains no directed edges or self loops 
(edges connecting a vertex to itself). It will be useful to represent the connectivity of the graph with 
the adjacency matrix, A = (Aij), where Aij = 1 if and only if ji G E (this "backwards" definition is 
not s t andar d, but it is useful for the linear algebra to come). The unfamiliar reader is directed to 



West! (|200ir ) for a thorough introduction to graph theory. 

When we wish to analyze the structure of a given network, one of the first things to examine 
is its degree distribution. The degree ki of a vertex i £ V is the number of edges incident to it, 
regardless of direction. The sequence fci, k2, ■ . . , kjy is called the degree sequence. In simple graphs, 
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this is just the size of a vertex's neighborhood. In more complicated graphs we usually speak instead 
of a vertex i's in-degree k\ , out-degree kf° , and undirected degree , all defined in the obvious 



way (jDodds et al 



201l[ ) . The degree distribution pk is a probability distribution which tells us how 
edges are distributed among nodes. The average degree fc avg = X/So ^Pk characterizes the overall 
density of edges. A degree distribution which is sharply peaked about its mean indicates a relatively 
homogeneous network where vertices tend to have the same number of incident edges. In contrast, 
a skewed distribution could result in some nodes having very high degree while the majority have 
low degree. This is common in real networks, and the tail of the degree d istribution often follows 



an approximate power law pk ~ k a for some exponent a (| Newman 



20031 ). 



2.2 Random graphs 

Random graphs are a family of models used to represent networks from the real world, although 
their suitability as such is questionable for reasons that will be elaborated below. Nevertheless, they 
are well-suited to analysis and as null models. 



The simplest random graph is the binomial model introduced by Erdos and Renyi (jBollobas 



200l|). Define a probability space Q{N,p) of graphs on N vertices where any of the (^) possible 



edges are chosen independently with probability p. The expected average degree in the network will 
be fc aV g = p(N— 1). We will often make statements about this network in the limit of large system size 
N — > oo. (In physics, this is referred to as the thermodynamic limit. Many expressions are simplified 
in performing this approximation, and derived results can be quite accurate given a large enough, 
albeit finite, N.) The degree distribution is a binomial distribution pk = ( N ~ 1 ^p k (l — p) N ~ 1 ~ k 
which is well- approximated by the Poisson distribution pk ~ (fc aV g) fc exp(— fc avg )/fc! with parameter 
fcavg = Np when TV is large and fc avg fixed. For this reason, Q(N,p) is also called the Poisson random 
graph model. As noted previously, many networks observed in the real world have heavy-tailed 
distributions, so the Poisson model is not suitable for those types of networks. Furthermore, real 
networks often contain a high density of triangles or clustering — friendship tends to be transitive. 
Poisson random graphs, on the other hand, have zero clustering in the thermodynamic limit. 

A mor e flexible generali z ation is the configuration model. In this mo del, either the degree 



sequence (|Mollov and Reed 



1995 



is given in advance (also see 



Newman 



19981) or t he expected degree sequence (jChung and Lu 



200J) 



20031 ); the models produce similar networks with subtle 



differences. The configuration model of can be thought of as a random wiring process as follows. 
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First, draw a degree sequence k\, . . . , fcjv independently according to the desired degree distribution 
Pk- We assign ki stubs (half-edges) to each vertex i € V. Next, choose a pair of stubs at random, 
connect them, remove both from the queue of stubs, and continue until all stubs are connect ed. 
Finally, pairs of edges arc chosen and their endpoints shuffled in the manner of . 



Milo et al 



(2002J) to 



ensure uniform sampling from all graphs with the given degree sequence. The resulting graph will 
have the imposed degree sequence, so long as the sum of the degrees is even. There may be self-loops 
or repeated edges, neither of which are allowable for simple graphs. However, we expect these to 
occur increasingly rarely for large TV, and they can be removed without affecting the re sulting graph 



much. Deg ree-degree correlat ions can be introduced by shuffling edges as described in 



Melnik et al 



2011 


) and 


Pavne et al. 


(2011) 



(|201lf ). The configuration model also lacks triangles and higher-order cliques 



in the thermo dynamic l i mit, a lthough some work has been done to create random graph models with 



clustering (see 



Hackett 



2011 



§4-5). 



2.3 Dynamical processes on networks 

Networks provide a structure on top of which all kinds of dynamical processes may take place. 
In many cases, the structure of the network itself heavily influences the dynamics. For instance, 
ferromagnetic materials can be modeled as atoms in a lattice network, where the bulk magnetization 
is the result of interactions between individual atoms' magnetic moments; strikingly different results 
occur in different dimensional lattices. Another example is a food web, a network of species connected 
by trophic interactions (who eats whom). The populations of the species in the ecosystem can be 
described by dynamical equations that reflect the structure of the network, and this can affect 
ecosystem stability. These are just two examples where dynamical processes on networks are a 
reasonable way to model sys tem behavior. Sin ce networks are ubiquitous structures, these models 



appear across all disciplines (jVespignani 



20121 ) 



2.3.1 Random Boolean networks 

Consider Boolean or binary state dynamics, where each node can be either "on" or "off" (in various 
contexts this can mean active/inactive, infected/susceptible, or spin up/spin down). The state of 
node i is encoded by a variable Xi G {0, 1}, and the system state is x = (xt). At each time step, nodes 
receive input from their neighbors in the (undirected) network. They then compute a function of that 
input, i.e., fi(x h , Xj 2 , . . . , Xj k . ) where ji, ... are the neighbors of node i. This determines their 



state in the next time step. Because the state space {0, 1}^ is finite, all trajectories are eventually 
periodic. The detailed structure of the cycles depends on the specific details of the network and the 
Boolean functions. 

The system described above, for ge neral fj : {0, l| fc ^ — > {0, 1}, is known as a Boolean network. 



Boolean networks were first studied by 



cells. See the review by 



Aldana et al 



Kauf fmanl (|1969[ ) as a model for dynamical behavior within 



20031 ). Most researchers have considered the dynamics on 



random i^-regular graphs with a parameter p that determines the bias between 0s and Is in the 
output of the update functions fi, which arc otherwise randomly chosen Boolean functions. 

As mentioned, deterministic Boolean network models must be eventually periodic. However, 
the behavior of the transient and the structure of the basins of attraction are different for different 
parameters K and p. In particular, ther e is a critical value o f the connectivity K c (p) that separates 



the transient dynamics into two phases (jAldana et al 



20031 ): 



1. Frozen, K < K c : The distance between nearby trajectories x(i) and x'(t) decays exponentially 
with time. 

2. Critical, K = K c : The temporal evolution of distance between trajectories is determined by 
fluctuations. 

3. Chaotic, K > K c : The distance between nearby trajectories grows exponentially with time. 



2.3.2 Social models 

When modeling social systems with a Boolean network, the nodes represent people and their states 
encode whether or no t they partic i pate i n a behavior, possess a certain b elief, etc. This could be 



riotin g or not rioting ([Granovetter 



19781 ). buying a particular style of tie (jGranovetter and Soong 



1986|), liking a particular band or style of music, or believing in some unintuitive or controversial 
idea, e.g. climate change. The state can represent any behavior with only two mutually exclusive 
possibilities. 

The function fi t hat determines how n ode i changes state is called its response function in 



sociological contexts. 



Schellind (11971 



19731 ) and 



Granovetterl (|l978l ) pioneered the use of threshold 



response functions in models of collective social behavior (although they were not the first; see 
the citations in their papers). This was based on the intuition that, for a person to adopt some 
new behavior, the fraction of the population exhibiting that behavior might need to exceed some 



G 



critical value, the person's threshold. (Mathematically, a threshold response function f(cj>; </> on ) with 



threshold 



returns if 6 < 



and 1 if <b > 



where the dynamics take place on a network by 
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Thes e models were generalized to the case 



Watta (|2002[ ). If we initialize a social network with 



some fraction of active nodes, some of their neighbors' thresholds may be exceeded and the activity 
can spread (depending on the distribution of thresholds and network structure). 

Standard threshold models have simple dynamical behavior: in a word, "spreading." If a single 
activation, on average, leads to more than one subsequent activation, then the spreading will be 
successful. The activity will increase in a sigmoid fashion until some final fraction of the network 
is activated. Th ese kinds of spreading are often studied using branching processes. See the book 

by 



Harris! (119631 ) for an overview of branching processes and the widely-used generating function 
formalism. Branching processes have been used to model, among other things, extinction of families, 
species, and genes, neutron cascades (as happens during nuclear chain reactions), and high energy 
particle shower s caused by cosmic rays. Th e generating function formalism developed in part by 



Newman (as in 



Newman 



2003 



Watts 



2002) to analyze spreading on networks is a straightforward 
application of the classical theory of branching processes. 

We note that threshol d random Boolean network models have been studied for the purpose of 



modeling neural networks (jAldana et al 



20031 ). Those models take place on signed, weighted graphs, 
which differ from the networks considered here, and the problems considered arc different. It would 
be interesting to explore the connections between our on-off threshold model (Section [3]) and other 
threshold Boolean network models. 



2.4 Some notation 

The Bachmann-Landau asymptotic notations are used throughout this thesis. When used carefully, 
asymptotic notation greatly improves the readability of analytic statements and pr oofs. It is also 
widely used in probability. For an overview of the notation's history and us age, see 



and the references therein. The notations we have used here are [citing from 



Knuthl dl97J)]: 



Knuthl (119761) 



0{f{n)) is the set of all g{n) such that there exist positive constants C and no with \g(n)\ < 
Cf(n) for all n>riQ. 



1 The edge case could be defined differently, but this will not influence the dynamics except in carefully constructed 
"pathological" scenarios. 
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• f2(/(n)) is the set of all g(n) such that there exist positive constants C and no with g(n) > 
Cf(n) for all ra > no- 

• o(/(n)) is the set of all g(n) such that g(n)/ f(ri) —> as n — > oo. 

• <?(n) ~ if g(n)/ f(n) — > 1 as n — >• oo. 

Formally, each of the above define sets of functions, but we often use statements such as "/ is O(g)" 
(read as "/ is big-oh of g") or "/ = 0(g)" to mean "/ £ 0(g)." 

When expressing probabilities we will use notation of the form P(x), P(x\y), etc. Here, P(x) is 
the probability that the random variable associated with x equals the specific value x. We leave out 
the random variables to avoid introducing unnecessary clutter. 

Symbols in boldface represent vector quantities or vector-valued functions, e.g., x = (xi). Sub- 
scripts have been left out in places for clarity. 



3 The on-off threshold model 



Here we study a s i mple extension of the classical thresh old models (such as 



Granovetter , 


1978 


Watts. 


2002 


Dodds and Watts. 


2004 



Schellind . 



1971 



1973 



20041 . among others): the response function also 



includes an off-thres hold. See Figure [2] for an exa mple on-off threshold response function. This is 



exactly the model of 
the following (also see 



Granovetter and Soong 



Granovetter and Soone ( 1986 ). but on a network. We motivate this choice with 



19861 ). (1) Imitation: the on state becomes favored as 



the fraction of active neighbors surpasses the on-threshold (bandwagon effect). (2) Non-conformity: 
the on state is eventually less favorable with the fraction of active neighbors past the off-threshold 
(reverse bandwagon, snob effect). (3) Simplicity: in the absence of any raw data of "actual" response 
functions, which are surely highly context-dependent and variable, we choose arguably the simplest 
deterministic functions which capture imitation and non-conformity. 

Let Sf = (V, E) be a graph with N = \ V\. Assign each vertex i £ V an on-threshold </> on ,i and an 
off-threshold <^ ff,j with < (f> n,i < 4>oS,i < 1- Then that node's response function fi(4>i', 0on,ij 4>oS,i) 
is 1 if O n,i < 4>i < '/'off,! and otherwise. Let x(0) £ {0, 1}^ be the initial states of all nodes. At 
time step t, each node i computes the fraction 4>i(t) of their neighbors in Sf who are active and takes 
the state Xi(t + 1) = /; (<j>i(t); 4>on,i, 4>oS,i) at the next time step. The above defines a deterministic 
dynamical system for a fixed graph and fixed thresholds. 
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We now make some quick remarks about the on-off threshold model. First, our model is a 
particular kind of Boolean network (Scction l2.3.1[) . Note that each node reacts only to the fraction 
of its neighbors who are active, rather than the absolute number, and the input varies from to 1 
in steps of 1/ki, where ki is node i's degree. Note that if 4> on ^ = the node activates spontaneously, 
and if </> ff.i = 1 we have the usual kind of threshold response function (without an off-threshold). 
A crucial difference between our model and many related thre shold models is that, in those 



Gleeson and Cahaland (|2007l ) 



models, an activated node can never reenter the susceptible state 
call this the permanently active property and elaborate on its importance to their analysis. Such 
models must eventually reach a steady state. When the dynamics are deterministic, this will be 
a fixed point, and in the presence of stochasticity the steady state is characterized by some fixed 
fraction of active nodes subject to fluctuations. The introduction of the off-threshold builds in a 
mechanism for node deactivation. Because nodes can now recurrently transition between on and off 
states, the deterministic dynamics can exhibit a chaotic transient (see Section l2.3.ip . and the long 
time behavior can be periodic with potentially high period. With stochasticity, the dynamics can 
be truly chaotic and never repeat. 

In the rest of this Section, we will describe how this model is different from the random Boolean 
networks in the literature. This is mainly due to the on-off threshold response functions we consider, 
but also the type of random graph on which the dynamics take place, varying amounts of stochasticity 
which we introduce in the networks and response functions, and the possibility of asynchronous 
updates. 



3.1 The networks considered 

The mean field analysis in Section [5] is applicable to any network which can be characterized by 
its degree distribution. As mentioned before, the vast majority of the theory of random Boolean 
networks assumes a regular random graph. Fortunately, such theories are easily generalized to 
other types of graphs with independent edges, such as Poisson and configuration model random 
graphs. Some specific results are given for the Poisson random graph Q(N, fc avg /iV), and these are 
the networks considered in Section [6] 
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Figure 2: An example on-off threshold response function. Here, (f> on = 0.33 and </> Q ff = 0.85. The 
node will activate if (j> on < (f> < 4> g, where </> is the fraction of its neighbors who are active. Otherwise 
it turns off. 



3.2 Stochastic variants 

The specific graph and node thresholds determine exactly which behaviors are possible. These are 
chosen from some distribution of graphs, such as Q(N, fc avg /iV), and some distribution of thresholds, 
given by the j oint density -P ( t6 n n, <? W). The specific graph and thresholds define a realization of 



the model (see 



Aldana et al 



2003J). When these are fixed for all time, we have, in principle, full 
knowledge of the possible model dynamics. Given an initial condition x(0), the dynamics x(t) are 
deterministic and known for all t > 0. 

With the introduction of noise, the system is no longer eventually periodic. Fluctuations at the 
node level allow a greater exploration of state space, and the behavior is comparable to that of the 
general class of discrete-time maps. Roughly speaking, the mean field theory we develop in Section[5] 
becomes more accurate as we introduce more stochasticity. 

We introduce randomness in two parts of the model: the network and/or the response functions. 
Allowing for the network and responses to be either fixed for all time or resampled each time step 
and taking all possible combinations yields four different designs (see Table [T]). 
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Rewiring network 


Fixed network 


Probabilistic response 


P-R 


P-F 


Deterministic response 


D-R 


D-F 



Table 1: The four different ways the model can be realized. These are the combinations of fixed or 
rewired networks and probabilistic or deterministic response functions. In the thermodynamic limit 
of the fully stochastic version (P-R), where the graph and response functions change every time step, 
the mean field theory is exact (see SecEJ. 

3.2.1 Rewired graphs 

First, the network itself can change every time step. This is the rewiring (R), as opposed to fixed 
(F), network case. For example, we could draw a new graph from G(N,k aV g/N) every time step. 
This amounts to rewiring the links while keeping the degree distribution fixed, and it is alternately 
known as a mean fi eld, "annealed" , or random mixing variant of the fixed or "quenched" model 



(Aldana et al 



2003J) 



3.2.2 Probabilistic responses 

Second, the response functions can change every time step. This is the probabilistic (P), as opposed 
to the deterministic (D), response function case. Again, we will need a well-defined distribution 
P(0 O n,0off) f° r the thresholds. This amounts to having a single response function, the expected 
response function 

f{4>)= /#on J #off P{4>oil, 0off)/(0; 0on, </>off)- (3.1) 

We call / : [0, 1] — > [0, 1] the probabilistic response function. Its interpretation is the following. 
For an updating node with a fraction <fi of active neighbors at the current time step, then, at the 
next time step, the node assumes the state 1 with probability f{4>) and the state with probability 

1-/(0)- 

3.2.3 The concept of "temperature" in the system 

In this thesis, the network and response functions are either fixed for all time or resampled every 
time step. One could tune smoothly between the two extremes by introducing rates at which these 
reconfigurations occur. These rates are inversely related to quantities that behave like temperature 
(one for the network and another for the response functions). For the fixed case, the temperature 
is zero, and there are no fluctuations, while in the stochastic case, the temperature is very large or 
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infinite, and fluctuations occur every time step. 



3.3 Synchronicity of the update 

Finally, we introduce a parameter a for the probability that a given node updates. When a — 1, all 
nodes update every time step, and the update rule is said to be synchronous. When a = 1/N, only 
one node is expected to update with each time step, and the update rule is said to be effectively 
asynchronous. This is equivalent to a randomly ordered sequential update. For intermediate values, 
a is the expected fraction of nodes which update each time step. 

4 Fixed networks 

Take the case where the response functions and graph are fixed (D-F), but the update may be 
synchronous or asynchronous. Let Xi(t) be the probability that node i is in state 1 at time t. and 
let fi{4>) = fi{4>',4>on,i,4>oS,i)- The dynamics follow the master equation 



Xi(t + 1) = afr [ V ] + (1 - a) Xi (t), (4.1) 



which can be written in matrix-vector notation as 

x(t + 1) = of (Tx(t)) + (1 - a)x(t). (4.2) 

Here T = D~ 1 A is sometimes called the transition probability matrix (in the context of a random 
walker), D is the diagonal degree matrix, and f = Note that if a = 1 we recover the fully 

deterministic response function dynamics, and Xi(t) = or 1 for all t. 

4.1 Asynchronous limit 

Here, we show that when a w 1/-/V, time is effectively continuous and the dynamic s can b e described 



by an ordinary differential equation. This is similar to the analysis of 



Gleeson 



(j2008f) . Consider 



2 Eqns. I I4.1H and 114.211 are not entirely correct when there are isolated nodes. In that case, ki = Ay = for 
certain i, thus the denominator in 114.111 is zero and D^ 1 undefined. If the initial network contains isolated nodes, we 
set all entries in the corresponding rows of T to zero. 
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Eqn. 14.21 Subtracting x(t) from both sides and setting Ax(£) = x(i + 1) — x(t) and Ai = 1 yields 

Ax(i) 



At 



= a(f(Tx(t))-x(t)). (4.3) 



Since a is assumed small, the right hand side is small, and thus Ax(t) is also small. Making the 
continuum approximation dx(t) jdt m Ax(i) / At yields the differential equation 

rl-x 

— = a(f(Tx)-x). (4.4) 

The parameter a sets the time scale for the system. From their form, similar asynchronous, contin- 
uous time limits apply to the dynamical equations in the densely connected case, Eqn. (|4.5[) . and in 
the mean field theory, Eqns. (|5.3[) and ([57 



4.2 Dense network limit for Poisson random graphs 

The following result is particular to Poisson random graphs, but similar results are possible for 
other random graphs with dense limit s. The norm alized Laplacian matrix is defi ned as C = I 



D-^AD- 1 / 2 , where lis the identity JWesd, [200lJ) . So T = D~ l / 2 {I - QD 1 / 2 . By 



Ohveiral d201dh . 



when fc avg is f2(logiV) there exists a typical Laplacian matrix C typ = In — ljsrl N /N [we let Ijv 
denote the length- N vector of ones and (•) the matrix transpose] such that the actual C « C typ 
in the induced 2-norm (spectral norm) with high probability. In this limit, if we assume uniform 
degrees fc, = k avg for all i £ V, then T w T typ = l?fV N /N. So T effectively averages the node states: 
Tx(£) w T typ x(i) = Yli=i X i(t)/N = 4>(t). Without a subscript, 4>{t) denotes the active fraction of 
the network at time t. We make the above approximation in Eqn. 14.21 and average that equation 
over all nodes, finding 



<f>{t + 1) = afm)) + (1 - a)<t>(t) = §(0(t); a, /), (4.5) 

where we have assumed that N is large and the average of nodes' individual response functions 
^2^ = ifi/N converges in a suitable sense to the stochastic response function /, Eqn. (|3.ip . This 
amounts to assuming a law of large numbers for the response functions, i.e., that the sample average 
converges to the expected function. Note that a tunes between the probabilistic response function 
$(</>; 1) = f(4>) and the 45° line $(0; 0) = <fi. Also, the fixed points of $ are fixed points of /, but 
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their stability will depend on a. 

When the network is dense, it ceases to affect the dynamics, since each node sees a large number 
of other n odes. Thus the network i s effec tively the complete graph. In this way we recover the map 
models of iGranovetter and SoonsJ (|1986h . 



5 Mean field theory 

In physics, making a "mean field" calculation refers to replacing the complicated interactions among 
many particles by a single interaction with some effective external field. There are analogous tech- 
niques for understanding networks dynamics. Instead of considering the \E\ interactions among the 
N nodes, network mean field theories derive self-consistent expressions for the overall behavior of the 
network, after ave raging over large sets of nod es. These have been fruitful in the study of random 



non-random (Melnik et al 



Boolean networks (Derrida and Pomeau 



19861 ) and can be surprisingly effective when networks are 

20111 ). 

We derive a mean field theory, in the thermodynamic limit, for the dynamics of the on-off 
threshold model by blocking nodes according to their degree class. This is equivalent to nodes 
retaining their degree but rewiring edges every time step. The model is then part of the well- 
known class of random mixing models with non-uniform contact rates. Probabilistic (P-R) and 
deterministic (D-R) response functions result in equivalent behavior for these random mixing models. 
The important state variables end up being the active density of stubs. In an undirected network 
without degree-degree correlations, the state is described by a single variable p(t). In the presence 
of correlations we must introduce more variables, i.e., Pk(t) , pv (t) , to deal with the relevant 
degree classes. 



5.1 Undirected networks 

To derive the mean field equations in the simplest case — undirected, uncorrelated random graphs 
— consider a degree k node at time t. The probability that the node is in the 1 state at time t + 1 
given a density p of active stubs is 



F k (p;f) = J2[.)p j a-p) k - j fU/k), 

3=0 Kh 



(5.1) 
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where each term in the sum counts the contributions from having 0, 1, . . . , k active neighbors. 
Now, the probability of choosing a r andom stub whic h ends at a degree k node is q k = fcpfe/fc avg 



in an uncorrelated random network ( Newman 



2003j). This is sometimes called the edge-degree 



distribution. So if all of the nodes update synchronously, the active density of stubs at t + 1 will be 



g(p;p k , f) = J2 qkF k (p; f) = £ ^F k (p; /). 
k=i fc=i Kav s 



(5.2) 



Finally, if each node only updates with probability a, we have the following map for the density of 
active stubs: 

p(t + l) = ag(p(t); Pk ,f) + (l-a)p(t) = G(p(t); Pk ,f,a). (5.3) 
By a similar argument, the active density of nodes is given by 



ct>{t + l)=a h(p(t); Pk , /) + (!- a)<P{t) = H{p(t),<l>(t);p k ,f, a), 



(5.4) 



where 



h(P:Pk, f) = ^2p k F k (p; /). 



fc=0 



(5.5) 



Note that the edge-oriented state variable p contains all of the dynamically important information, 
rather than the vcrtcx-oricntcd variable <b. 



5.2 Analysis of the map equation 



The function F k (p; f) is known in 



(in the variable p) of / (jPhillips 



; )olynomial approx imation theory as the fcth Bernstein polynomial 



2003 



Pena 



19991 ). These are approximating polynomials which 



have applications in computer graphics due to their "shape-preserving properties." The Bernstein 
operator takes / h- )• F k . This is a linear, positive operator which preserves convexity for all k and 
exactly interpolates the endpoints /(0) and /(l). Immediate consequences include that each F k is a 
smooth function and the kth derivatives Fj^\x) — > f^(x) where f^ k \x) exists. For concave / (such 
as the tent or logistic maps), we have concave F k for all k and F k /* f uniformly. This convergence 
is typically slow. Importantly, F k /* f implies that g(p;p k , /) < / for any degree distribution p k . 

In some cases, the dynamics of the undirected mean field theory given by p(t + 1) = G(p(t)), 
Eqn. (|5.3|) . are effectively those of the map from the dense limit Eqn. (|4.5p . We see that g, 
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Eqn. (|5.2p . can be seen as the expectation of a sequence of random functions Fk under the edge- 
degree distribution qk (indeed, this is how it was derived). From the convergence of the FJt's, we 
expect that g(p\Pk,f) ~ f(p) if the average degree fc avg is "large enough" and the edge-degree 
distribution has a "sharp enough" peak about fc avg (we will clarify this soon). Then as fc avg — > oo, 
the mean field coincides with the dense network limit we found for Poisson random graphs, Eqn. 14.51 
Some thought leads to a sufficient condition for this kind of convergence: the standard deviation 
c (A: aV g) of the degree distribution must be o(/c avg ). In Appendix lAl we prove this as Lemma [1] 

In general, if the original degree distribution pk is characterized by having mean fc avg , variance 
a , and skewness 71, then the edge-degree distribution qk will have m.6ctll Aj aV g "h^^/^avg and Variance 
<r 2 [l +7icr/£; aV g — (er/£; aV g) 2 ]. Considering the behavior as k avg oo, we can conclude that requiring 

a G o(fcavg) and 71 G o(l) are sufficient conditions on pk to apply Lemma [TJ Poisson degree 

1/2 

distributions (a = -^/fcavg and 71 = fc avg ) fit these criteria. 

5.3 Generalized random networks 

In the most general kind of random networks, edges can be undirected or directed, and we then 
denote node degree by a vector k= (fc( u \ fc«, fc^) t . The de grce distribution is written aspk = -P(k). 
There may also be correlations between node degrees. Correlations of this type are encoded by the 
conditional probabilities 

P]lIl' = Pfe, undirected I k') 
Pkk/ = P(k, incoming|k') 
Pk°k' = -P(k,outgoing|k'), 



the probability that an edge starting at a degree k' node ends at a degree k node and is, respectively, 
undirected, incoming, or out going relative to the destination degree k node. We introduced this con- 



vention in a series of papers (jPavne et al 



2011 



Doddsetal 



20111 ) . These conditional probabilities 



can also be defined in terms of the joint distributions of node types connected by undirected and 
directed edges. The mean field equations for the on-o ff threshold model are closely related to the 

Eqns. (13-15)). We omit 



equations for the time evolution of a contagion process (jPavne et al 



2011 



a detailed derivati on, since it is similar to that in Section [5. II fsee also 



Gleeson and Cahalane 



2007 



Payne et al 



20111 ) . The result is a coupled system of equations for the density of active stubs which 
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now may depend on node type (k) and edge type (undirected or directed): 



(u) 



(*+i)=«£*& EE 



*™ ^ /fcH'x / fc (i)' 



3*=0ji=0 



3u J \ 3 



ji r 



+ (l-a)/>L u) (t) 



i^'-n) ( j u +ji 



k' j„=0j»=0 \ 3 U / \ J 1 



+ (i-«)Pk°W. 



fcW + fcW 



(5.6) 



(5.7) 



The active fraction of nodes at a given time is given by: 



^(t+i) =«E^ k E E 



k j„=Oj«=0 



Ju 



+ (1 - a)0(t) . 



(fc (,) -J.) 

/ 



J« + 3i 

jfe(u) + jfe(i) 



(5. 



6 Poisson random graphs with tent map average response 
function 

The results so far have been entirely general, in the sense that the underlying network and thresholds 
are arbitrary. Now we apply the general theory to the case of Poisson random graphs with a simple 
distribution of thresholds. 

The networks we consider are Poisson random graphs from G(N,k avg /N). The thresholds <p on 
and o ff arc now distributed uniformly on [0, 1/2) and [1/2, 1), respectively. This distribution results 
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Figure 3: The tent map probabilistic response function f(p), Eqn. (|6.1|) . Note that we use the 
argument p for comparison with the edge maps g(p; fcavg) = g(p;pk, /), Eqn. (|5.2j) . shown for fc avg = 
1, 10, 100. These pk are Poisson distributions with mean k avs . As fc avg increases, g(p; k avg ) increases 
to f(p). 



in the probabilistic response function (sec Figure [3]): 



f(4>) = 



20 if O<0<l/2, 
2-20 if 1/2<0<1. 



(6.1) 



19961 ). We thus expect 



The tent map is a well-known chaotic map of the unit interval (jAlligood et al. 
the on-off threshold model with this probabilistic response function to exhibit similarly interesting 
behavior. 



6.1 Analysis of the dense limit 

When the network is in the dense limit (Section 14. 2[) . the dynamics follow <f>(t + 1) = <I>(0(i); a), 
where 

{l + a)4> if O<0<l/2, 

(l-3a)0 + 2a if 1/2 < < 1. 



$(0;a) = a/(0) + (l-a) 



(6.2) 



Solving for the fixed points of ■£(</>; a), we find one at = and another at = 2/3. When a < 2/3, 
the nonzero fixed point is attracting for all initial conditions except = 0. When a = 2/3, [1/2, 5/6] 
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■0.1 



update probability a 

Figure 4: Bifurcation diagram for the dense map ( &(</'; a), Eqn. (|6.2[) . This was generated by iterating 
the map at 1000 a values between and 1. The iteration was carried out with 3 random initial 
conditions for 10000 time steps each, discarding the first 1000. The 0-axis contains 1000 bins and 
the invariant density, shown by the grayscale value, is normalized by the maximum for each a. With 
a < 2/3 (not shown), all trajectories go to the fixed point at </> = 2/3. 

is an interval of period 2 centers. Any orbit will eventually land on one of these period 2 orbits. 
When a > 2/3, this interval of period 2 orbits ceases to exist, and more complicated behavior ensues. 
Figure HI shows the bifurcation diagram for $(0; a). From the bifurcation diagram, the orbit appears 
to cover dense subsets of the uni t interval when a > 2 /3. The bifurcation diagram appears like that 



of the tent map (not shown; see 



Alligood et al 



1996 



Doddsetal 



20121 ). except that the branches 



to the right of the first bifurcation point are separated here by the interval of period 2 orbits. 



6.1.1 The effect of conformists, an aside 

Suppose some fraction c of the population is made up of individuals without any off-threshold 
(alternatively, each of their off-thresholds <f> s = 1). These individuals are conformist or "purely 
pro-social" in the sense that they are perfectly happy being part of the majority. For simplicity, 
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assume a = 1. The map $(</>; c) = 20 for < <p < 1/2 and 2 - 2(1 - c)0 for 1/2 < (j> < 1. If c> 1/2, 
then the equilibrium at 2/3 is stable. Pure conformists, then, can have a stabilizing effect on the 
process. We expect a similar effect when the network is not dense. 



6.2 Mean field 



6.2.1 Analysis 

In this specific example, w e can write t he degree-dependent map Ff. (p; /) in terms of incomplete 



regularized beta functions (|NIST 



20121 ). Since / is understood to be the tent map, we will write 



Fk{p\ f) = Fk{p)- First, use the piecewise form of Eqn. (|6.1[) to write 



M 



3=0 



j=M+l 



= 2 -2p- 2^ (1 - P f-> + (i) £ (1 - 

We have let M = \k/2\ for clarity ([-J and [•] are the floor and ceiling functions) and used the fact 
that the binomial distribution (^)pi(l — p) fc_ -' sums to one and has mean kp. For n < M, we have 
the identity 



E C?')« • VX 1 - ^ = P n (fc)™/i-p(fc - M, A/ - n + 1) 

3=0 



(6.3) 



where I x {o-, b) is the regularized incomplete beta fu nction and (fc)„ = fc(fc — 1) • • • (k — (n — 1)) is the 



falling factorial (IWinkler et al 



1972 



NIST 



20121 ). This is an expression for the partial (up to M) 



nth factorial moment of the binomial distribution with parameters k and p. Note that when n = 
we recover the well-known expression for the binomial cumulative distribution function. So, 



F k (p) =2p - 4p7 p (M, k - M) + 2I p (M + l,k- M) (6.4a) 
= (2 - 2p) - (27!_ p (fe - M, M + l)- Aph- P {k - M, M)). (6.4b) 

When < p < 1/2, the form of Eqn. ()6.4a|) can be used to show directly that Fk(p) is bounded 
above by the tent map f{p), which we already knew from the properties of the Bernstein polynomials. 
A similar approach works for the region 1/2 < p < 1 using Eqn. (|6.4b|) . We find a weak bound for 
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the rate at which Fk(p) converges to f(p). For p < 1/2, using Eqn. (|6.4a[) . 



f(p) - F k {p) = 4pI p (M, k - M) - 2I P (M + 1, k - M) 

< 2 (I P (M, k - M) - I p (M + l,k- M)) 
p M (l-p) k - M 



< 2 



< 



MB(M, k — M) 
2(l/2) fc 



MB(M, k - M) 

4r(fc) 

^ WmW ( even } 

i r(fe/2 + i/2) 1/2 
- * r(*/2) ~ U{ h 



NIST1 I2OIJ). We find the same Oik" 1 / 2 ) behavior using 



where we have used identities from 
Eqn. (|6.4bp in the region p > 1/2. 

Finally, note that the active edge fraction p(t) w </>(£), the active node fraction. This is because 
Poisson random graphs are highly regular, with = kpk/k avg = pk-i ~ Pk- Thus the mean field 
dynamics for active edge density are effectively the same as for active node density. 

6.2.2 Numerical algorithm 

The map g(p; pk, f) is parametrized here by the network parameter fc a vg, since pk is fixed as a Poisson 
distribution with mean /c avg and / is the tent map, and we write it as simply g(p; fc avg ). To evaluate 
<?(/o;fc a vg), we compute Fk(p) using Eqn. (|6.4[) and constrain the sum in Eqn. (|5.2[) to values of k 
with [fcavg — 3 -\/fcavgJ < k < [fc avg + 3 y^fcavg] • This computes contributions to within three standard 



deviations of the average degree in the graph, requiring only 0(y/k avs ) evaluations of Eqn. 
The representation in Eqn. (|6.4[) allows for quick numerical evaluation of Fk(p) for any fc, which we 
performed in MATLAB using the built-in routines for the incomplete beta function. 

In Figure [31 we show g(p;k avg ) for /c avg = 1,10, and 100. We confirm the conclusions of Sec- 
tion [521 s(p;^avg) is bounded above by f(p), and g{p',k avg ) /*■ f(p) as fc avg — > 00. Convergence is 
slowest at p = 1/2, and the kink that the tent map has there has been smoothed out by the effect 
of the Bernstein operator. 
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6.3 Simulations 

We performed direct simulations of the on-off threshold model for the D-F, P-F, and P-R designs, in 
the abbreviations of Table [TJ Unless otherwise noted, N = 10 4 . For all of the bifurcation diagrams, 
the first 3000 time steps were considered transient and discarded, and the invariant density of p was 
calculated from the following 1000 points. For plotting purposes, the invariant density was normal- 
ized by its maximum at those parameters. For example, in Figure HI we plot P {4>\ot) / max^ P(<t>\a) 
rather than the raw density P(<p\a). 

To compare the mean field theory to those simulations, we numerically iterated the edge map 
p(t + 1) = G(p(t); fc avg ,a) for different values of a and k avg . We then created bifurcation diagrams 
of the possible behavior in the mean field as was done for the simulations. 



6.4 Results 

We show results for the deterministic model on a small network in Figure [5j Here, N = 100 and 
fcavg = 17. Starting from an initial active node at t = 0, the active population grows monotonically 
over the next 6 time steps. From t = 6 to t 80, the transient time, the active population fluctuates 
in a similar manner to the stochastic case. After t « 80, the state collapses into a period 4 orbit. 
We call the overall period of the system its "macropcriod." Individual nodes may exhibit different 
"microperiods." Note that the macroperiod is the lowest common multiple of the individual nodes' 
microperiods. In Figure [SJ we observe microperiods 1, 2, and 4 in the timcseries of individual node 
activity. A majority of the nodes end up frozen in the on or off state, with approximately 20% of the 
nodes exhibiting cyclical behavior after collapse. The focus of this thesis has been the analysis of 
the on-off threshold model, and the D-F case has not been as amenable to analysis as t he stochastic 
cases. A deeper examination through simulation of the deterministic case will appear in 



Doddsetal 



(|2012h 



We explore the mean field dynamics by examining the limiting behavior of the active edge fraction 
p under the map G(p; fc avg , a). We simulated the map dynamics for a mesh of points in the (fc avg , a) 
plane. We plotted the 3-dimensional (3-d; A-d denotes A-dimensional) bifurcation structure of the 
mean field theory in Figure [51 We also made 2-d bifurcation plots for fixed fc avg and a slices through 
this volume, shown in Figures [7] and [5] In all cases, the invariant density of p is normalized by its 
maximum for that (fc avg , a) pair and indicated by the grayscale value. 
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Figure 5: Deterministic (D-F) dynamics on a small graph. Here, N = 100 and fc avg = 17. On the 
left, we plot the state evolution over time. The upper plot shows individual node states (black = 
active) sorted by their eventual level of activity, and the lower plot shows the total number of active 
nodes. We see that the contagion takes off, followed by a transient period of unstable behavior until 
around time step 80, when the system enters a macroperiod 4 orbit. Note that individual nodes 
exhibit different microperiods (sec Sec. 16.41) . On the right, we show the network itself with the initial 
seed node in black in the lower right. 
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Figure 6: The 3-dimensional bifurcation diagram computed from the mean field theory. The axes X 
= average degree fc avg , Y = update probability a, and Z = active edge fraction p. The discontinuities 
of the surface are due to the limited resolution of our simulations. See Figure [7] for the parameters 
used. This was plotted in Paraview. 
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Figure 7: Mean field theory bifurcation diagram slices for various fixed values of fc avg and a. The top 
row (a-c) shows slices for fixed fc avg . As fc avg — > oo, the fc avg -slice bifurcation diagram asymptotically 
approaches the bifurcation diagram for the dense map, Figured Note that the first bifurcation point , 
near 2/3, grows steeper with increasing A: avg . The bottom row (d-f) shows slices for fixed a. The 
resolution of the simulations was a = 0.664, 0.665, . . . , 1, fc avg = 1, 1.33, . . . , 100, and p bins were 
made for 1000 points between and 1. 




Figure 8: Bifurcation diagram from fully stochastic (P-R) simulations. The same parameters were 
used as in Figure [71 which has the same structure only more blurred. 
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The mean field map dynamics exhibit period-doubling bifurcations in both parameters fc avg and 
a. Visualizing the bifurcation structure in 3-d (Figure [B]) shows interlacing period-doubling cascades 
in the two parameter dimensions. These bifurcations are more clearly resolved when we take slices 
of the volume for fixed parameter values. The mean field theory (Figure [7]) closely matches the 
P-R simulations (FigureH]). The first derivative dG(p; fc avg , a) /dp < d$(p; a) /dp for any finite fc avg , 
so the bifurcation point a = 2/3 which we found for the dense map $ is an upper bound for the 
first bifurcation point of G. The actual location of the first bifurcation point depends on fc avg , but 
a = 2/3 becomes more accurate for higher fc avg (it is an excellent approximation in Figures [7J; 
andJSfc, where fc avg = 100). When a = 1, the first bifurcation point occurs at fc avg ps 7. 

The bifurcation diagr am slices resemble each o ther and evidently fall into the same universality 



class as the logistic map (jFeigenbaum 



1978 



19791 ). This class contains all Td maps with a single, 
locally-quadratic maximum. Due to the properties of the Bernstein pol ynomial s , Fk(p \ f) will uni- 



2003J). So this will 



versally have such a quadratic maximum for any concave, continuous / (| Phillips , 
also be true for g(p; fc avg , /) with fc avg finite, and we see that fc avg partially determines the amplitude 
of that maximum in Figure [3] Thus fc avg acts as a bifurcation parameter. The parameter a tunes 
between G(p; fc avg , 1) = g(p; fc avg , /) and G(p; fc avg , 0) = p, so it has a similar effect. Note that the 
tent map / and the dense limit map $ are kinked at their maxima, so their bifurcation diagrams are 
qualitatively different from those of the mean field. The network, by constraining the interactions 
among the population, causes the mean field behavior to fall into a different universality class than 
the response function map. 



7 Conclusions 

We constructed the on-off threshold model as a simple model for social contagion resulting from 
limited imitation. We see that including an aversion to total conformity results in more complicated, 
even chaotic dynamics. This model also allows us to study the effects of differing amounts of fixedness 
in the social network and individual response functions, and we developed a detailed mean field 
theory which is exact for random mixing versions of the model. Finally, we applied the theory to a 
specific case, where the network is Poisson and the response functions average to the tent map. 

The model exhibits rich mathematical behavior. The deterministic case, which we have barely 
touched on here, merits further study. In particular, we would like to characterize the distribution of 



2G 



periodic sinks, how the collapse time scales with system size, and how similar the transient dynamics 
are to the mean field dynamics. 

Furthermore, the model should be tested on realistic networks. These could include power law 
or small world rand om graphs, or real social networks gleaned from data. In a manner similar to 



Melnik et al 



(|201l[ ). one could evaluate the accuracy of the mean field theory for real networks. 
Finally, the ultimate validation of this model would emerge from a better understanding of social 
dynamics themselves. Characteriza tion of people's true response functions is t herefore critica l (som e 



work has gone in this direction; see 



Centola 



2010 



2011 



Romero et al 



2011 



Ugander et al 



20121 ) 



Comparison of model output to large data sets, such as observational data from social media or online 
experiments, is an area for further experimentation. This might lead to more complicated context- 
and history-dependent models. As we collect more data and refine experiments, the eventual goal 
of quantifiably predicting human behaviors, including fashions and trends, seems achievable. 
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Appendix A Proof of Lemma 1 

Lemma 1. For k > 1, let ff. be continuous real-valued functions on a compact domain X with 
fk f uniformly. Let pk be a probability mass function on Z + parametrized by its mean \x and with 
standard deviation cr(fi), assumed to be o(fi). Then, 

lim V p k f k = f. 

l— ►no \ ^ ' / 



fl — ^ OO 



\k=0 



Proof. Suppose < a < 1 and let K = \ji — /i a J . Then, 

oo K oo 

fc=0 fc=0 k=K+l 

Since fk — > f uniformly as k — > oo, for any e > we can choose [i large enough that 

\f k (x)-f(x)\<e (A.2) 
for all k > K and all x € X. Without loss of generality, assume that \fk\ < 1 for all k. Then, 

2 



f\ < 



The a I \x a term is a consequence of the Chebyshev inequality dBollobasl . l200lh applied to the first 
sum in (|A.1[) . Since a grows sublinearly in [i, this term vanishes for some < a < 1 when we take 
the limit /i — > oo. The e term comes from using (|A.2[) in the second sum in (|A.1[) . and it can be 
made arbitrarily small. □ 



Appendix B Online material 

To better explore the 3-d mean field bifurcation structure, we created movies of the the /c avg and a 

slices as the parameters are dialed. The videos are available at 

http : //www . uvm . edu/~kharris/ on-of f -threshold/bif urc_movies . zip 

Also, a VTK file with the 3-d bifurcation data, viewable in Paraview is at 

http : //www. uvm. edu/~kharris/on-of f-threshold/volume_normalized . vtk.zip| 

Videos of the individual-node dynamics for small networks in the D-F and P-F cases are shown 
for some parameters which produce interesting behavior. These are available at 
|http : //www . uvm . edu/~kharr is/on- of f -thr eshold/graphji ovies . zip| 

The D-F, P-F, and P-R cases were implemented in Python. The code is available at 
|http : //www . uvm . e du/~kharris/on-off -threshold/code . zip 
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